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Section  1 


INTRODUCTION 


Under  Air  Force  Weapons  Laboratory  (AFWL)  contract  AF29  (601 ) - 
5399  IBM  developed  an  automated  transient  circuit  analysis  digital  computer 
program,  PREDICT,  for  analyzing  the  responses  of  electronic  circuits  in 
transient  radiation  environments,  and  a  system  of  digital  computer  programs, 
TREAT,  for  processing  component  radiation  effects  data.  These  programs 
were  subsequently  maintained,  modified,  and  distributed  to  various  aero¬ 
space  companies  throughout  the  United  States.  After  extensive  application 
of  PREDICT,  a  number  of  modifications  and  improvements  were  suggested. 

Under  AFWL  contract  AF29  (601)-6489  IBM  developed  and  evaluated  the 
suggested  improvements,  and  recommended  that  a  new  circuit  analysis  pro¬ 
gram  incorporating  these  features  be  developed.  The  effort  under  AFWL 
contract  AF29  (601)-6852  consisted  primarily  of  formulating  and  coding  this 
hew  circuit  analysis  program,  which  was  called  SCEPTRE.  The  most  recent 
contract,  F29601-70-C-0038,  was  awarded  for  the  purpose  of  implementation 
of  a  series  of  desirable  improvements. 
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Section  II 


FORMULATION  AND  THEORY 


2.  1  INTRODUCTION 


Any  automatic  transient  analysis  program  is  designed  to  relieve  the 
user  of  the  necessity  of  writing  and  programming  the  differential  and  algebraic 
equations  that  describe  networks.  PREDICT  and  other  programs  already 
perform  this  basic  task,  but  the  degree  of  flexibility  permitted  the  user  varies 
widely  among  programs.  SCEPTRE,  written  as  a  successor  to  PREDICT, 
incorporates  many  improvements. 

This  section  presents  the  formulation  and  theory  that  serve  as  the  basis 
for  SCEPTRE.  Therefore,  the  discussion  is  mathematically  oriented.  Those 
features  of  SCEPTRE  that  are  not  mathematical  are  not  included  here.  For 
example,  the  extremely  useful  features  of  Rerun  and  Model  Storage  are  not 
described. 

2.2  GENERAL  SOLUTION 


SCEPTRE  consists  of  two  separate  formulations  that  combine  to  pro¬ 
duce  the  general  transient  solutions  of  a  given  network.  One  is  referred  to 
as  the  initial -conditions  program.  This  will  determine  the  network  voltages 
that  prevail  before  any  time-varying  forcing  function  is  applied.  This  program 
does  not  treat  time  as  an  independent  variable;  instead  it  holds  time  constant, 
and  iterates  on  selected  voltages.  The  output  of  the  initial-conditions  pro¬ 
gram  may  be  obtained  independently  or  it  may  be  automatically  used  as  the 
starting  point  for  the  transient  program.  Thus,  the  output  of  the  initial- 
conditions  program  effectively  supplies  the  initial  conditions  for  the  system 
of  differential  equations  that  are  solved  by  the  transient  program. 

The  second  program  is  called  the  transient  program.  This  program 
uses  time  as  an  independent  variable  and  solves  systems  of  differential  equa¬ 
tions  as  functions  of  time.  The  output  of  this  program  represents  the  tran¬ 
sient  response  of  a  given  network.  As  implied  above,  the  transient  program 
may  be  used  in  conjunction  with  the  initial-conditions  program,  or  it  may  be 
used  by  itself  if  the  initial  conditions  of  the  network  are  known. 
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The  general  solution  procedure  described  here  concerns  the  definition 
of  the  terms,  matrices  and  procedures  that  are  common  to  both  programs. 
Other  parts  of  this  volume  will  provide  the  detailed  explanations  and  deri¬ 
vations.  The  first  stop  in  either  program  is  the  construction  of  a  tree* 
(Ref.  1)  according  to  prescribed  rules,  which  differ  for  the  two  programs. 
This  permits  formation  of  a  B  matrix  that  effectively  expresses  link  volt¬ 
ages  in  terms  of  tree  branch  voltages  and  tree  branch  currents  in  terms  of 
link  currents.  Figure  1  shows  a  composite  B  matrix  that  contains  all  pos¬ 
sible  element  classifications  and  submatrfces.  This  matrix  is  derived  in 
appendix  I.  The  element  classifications  are  given  in  table  I. 
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Figure  1.  Composite  B  Matrix  in  SCEPTRE 


Reproduced  from 

belt  available  copy^jji^ 


*A  tree  is  defined  as  any  connected  network  subgraph  that  contains  all 
nodes  of  the  network  but  no  complete  loops.  All  circuit  elements  that  are 
members  of  the  tree  are  termed  tree  branches.  All  circuit  elements  ex¬ 
cluded  from  the  tree  are  termed  links.  A  "C"  tree  is  defined  in  this  report 
ns  one  in  which  tree  members  are  chosen  in  the  preference  order  E,  C,  T! 
and  L.  All  current  sources  (J)  must  be  excluded  from  the  "C”  tree.  Thor» 
fore,  these  sources  are  links.  A  cut  set  is  defined  as  that  group  ok  eiemenls 
that  would  isolate  two  groups  of  nodes  when  removed  from  a  network. 
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Table  1 


ELEMENT  CLASSIFICATIONS  IN  SCEPTRE 


Class 

Element 

1 

Capacitor  links 

2 

Resistor  links 

3 

Inductor  links 

4 

Capacitor  tree  branches 

5 

Resistor  tree  branches 

6 

Inductor  tree  iches 

7 

Voltage  sources 

8 

Independent  current  sources 

9 

Primary  current  sources  (dependent  on  voltage  across  terminals). 
This  class  will  appear  only  in  the  derivation  of  the  initial- 
conditions  program. 

0 

Secondary  current  sources  (dependent  on  other  current 
sources).  This  class  will  appear  only  in  the  derivation  of  the 
initial- conditions  program. 

Y 

Voltage  sources  that  are  dependent  on  resistor  voltages.  This 
class  will  appear  only  in  the  derivation  of  the  transient  program. 

X 

Current  sources  that  are  dependent  on  resistor  currents.  This 
class  will  appear  only  in  the  derivation  of  the  transient  program. 

The  following  matrices  ard  vectors  may  also  be  defined  as  a  result  of 
the  element  classification  ir.  table  I: 


•  ~  a  diagonal  matrix  composed  only  of  resistor  links 

*  G22  "  R22 

•  -  a  diagonal  matrix  composed  only  of  resistor  tree  branches 

•  Gre.  -  Rr  r  ~ 1 

55  55 

»  -  a  diagonal  matrix  composed  only  of  capacitor  links 

•  sn  Si'1 

•  ^44  "  a  diagonal  matrix  composed  only  of  capacitor  tree  branches 

•  s44  =  C44 

•  Lgg  -  a  matrix  composed  only  of  inductor  links  and  the  mutual 

inductance  between  inductor  links 

•  -  a  matrix  composed  only  of  the  mutual  inductance  between 

inductor  links  and  inductor  tree  branches 

•  Lg6  -  a  matrix  composed  only  of  inductor  tree  breaches  and  the 

mutual  inductance  between  inductor  tree  branches 


Ggg  -  a  diagonal  matrix  composed  only  of  the  voltage  derivatives  of 
primary  current  sources 

Ij,  Vj  -  vectors  composed  only  of  the  currents  or  voltages  asso¬ 
ciated  with  capacitor  links 
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•  Ig*  V 2  -  vectors  composed  only  of  the  currents  or  voltages  asso¬ 

ciated  with  resistor  links 

•  Ig,  Vg  -  vectors  composed  only  of  the  currents  or  voltages  asso¬ 

ciated  with  inductor  links 

■J  i 

•  I^,  Vj  -  vectors  composed  only  of  the  currents  or  voltages  asso¬ 

ciated  with  capacitor  tree  branches 

©  Ig,  Vg  -  vectors  composed  only  of  the  currents  or  voltages  asso¬ 
ciated  with  resistor  tree  branches 

•  Ig  Vg  -  vectors  composed  only  of  the  currents  or  voltages  asso- 

’  dated  with  inductor  tree  branches 

•  Jg,  Vg  -  vectors  composed  only  of  the  currents  or  voltages  asso¬ 

ciated  with  independent  current  sources 

•  Jg,  Vg  -  vectors  composed  only  of  the  currents  or  voltages  asso¬ 

ciated  with  primary  current  sources 

•  Jg,  V  -  vectors  composed  only  of  the  currents  or  voltages  asso¬ 

ciated  with  secondary  current  sources 

•  E?  -  a  vector  composed  only  of  voltage  sources.  , . 

2.3  TRANSIENT  SOLUTION 

The  state  variables  of  any  system  can  be  defined  as  the  minimum  set  : 
of  quantities  that  will  suffice  to  determine  all  other  quantities  in  the  system 
at  any  instant.  It  can  be  shown  that  the  knowledge  of  the  set  of  all  capacitor 
tree  branch  voltages  (V4)  and  inductor  tree  link  currents  (I3)  is  sufficient  to 
determine  all  other  element  currents  and  voltages,  and  furthermore,  that  this 
selection  of  state  variables  will  allow  network  formulation  in  terms  of  first- 
order  differential  equations.  The  starting  point  of  the  derivation  then  is  that 
quantities  V4  and  I3  are  known  and  the  list  of  unknowns  is  made  up  of  Vi>  V2. 
V3,  V5,  V6>  II,  \2f  l4»  I5,  l6>  17  and  Vs.  In  addition,  the  derivatives  of  the 
state  variables,  V4  and  I3,  must  be  obtained  in  preparation  for  the  numerical 
integration  routine  which  produces  the  updated  state  variables  that  are  valid 
at  the  next  time  increment. 

Some  oi  the  equations  needed  to  solve  the  unknown  quantities  may  be 
obtained  from  the  transient-solution  B  matrix  (figure  2).  The  B  matrix 
itself  arises  from  a  "C"  tree,  which  is  formed  by  an  E,  C,  R,  L  preference 
order.  Note  that  the  B  matrix  differs  from  the  composite  matrix  of 
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Figure  2.  Transient -Solution  B  Matrix  in  SCEPTRE 


Figure  1 1  in  that  some  submatrices  are  zero  valued*  and  that  no  distinction 
is  made  between  types  of  sources.  Since  link  voltages  can  be  written  in 
terms  of  tree  branch  voltages  directly  from  the  B  matrix,  the  following  equa¬ 
tions  may  be  written  in  matrix  form: 


-B14  V4 


-B24  V4 


B34  V4 


-B84  V4 


-  b1?  E? 


-  B25  V5 


'  B35  V5 


“  B85  V5 


"  B27  E7 

"  B36  V6  '  B37  E7 
"  B86  V6  '  B87  E7 


(1) 

(2) 

(3) 

(4) 


Since  t'ree  branch  currents  can  be  written  in  terms  of  link  currents,  there 
arises: 

*4  =  B14  *1  +  B24  l2  +  B34  h  +  B84  J8  ^ 


♦For  example,  B^  must  be  zero  since  the  ”C"  tree  preference  prohibits  the 
possibility  of  a  capacitor  link  closing  a  loop  that  contains  a  resistor  tree 
branch. 
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(6) 


X5  = 


B25Tl2  +  B35Tl3  +  B85Tj8 


!6  = 


B36  l3  +  B86  J8 


h  =  B17  !1  +  B27  l2  +  B37  *3  +  B87  J8 


(7) 

(8) 


where  the  superscript  T  is  used  to  indicate  the  transpose  of  a  matrix. 
Two  additional  equations  may  be  obtained  from  differentiating  Equations 
(1)  and  (7)  yielding: 


• 

V1  ■ 

“B14  ^4  -  B17  E7 

(9) 

%  ■ 

B36T  ^3  +  B86  ^8 

(10) 

Note  that  source  derivatives  E7  and  jg  have  been  introduced  in  the  last  two 
equations.  *  These  equations  together  with  a  few  fundamental  relations  will 
be  used  to  derive  all  of  the  network  currents  and  voltages  in  terms  of  known 
quantities. 

2.  3.  1  SOLUTION  OF  RESISTIVE  QUANTITIES 

The  resistive  quantities  of  interest  are  I2,  I5,  V2 ,  and  Vg  .  The  funda¬ 
mental  voltage  current  relations  for  resistors  permit: 

V2  =  «22  h  (11) 

V5  '  R55  h  (12) 

12  may  be  explicitly  solved  for  by  manipulation  of  equations  (2),  (6),  (11)  and 
(12)  to  get 

^  -  V1  {  -®24  V4  '  B27  E7  -  B25  ^5  [B35T  !3  *  B85T  J8  ]}  <13> 

where 

MR  =  R22  +  B25  R55  B25 


♦See  section  2.  5  for  a  discussion  on  source  derivatives. 
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The  significance  of  equation  (13)  is  that  the  vector  of  resistor  link  currents 
has  been  solved  in  terms  of  all  known  quantities  since  the  right  side  of  the 
equation  is  composed  entirely  of  state  variables  V4  and  I3,  known  sources 
E7  and  Jg,  and  known  incidence  submatrices.  Once  lo  is  known,  the  vectors 
I5,  V«>,  and  V5  can  be  determined  from  Equations  (6),  (11).  and  (12), 
respectively. 

An  alternate  approach  may  sometimes  be  preferable.  Equations  (2), 
(6),  (11),  and  (12)  may  also  be  manipulated  to  obtain 

V5  *  V1  {-B25T  °22  [B24  V4  +  B27  E7]  +  B35T  J3  +  B85Tj8}  (15) 

where 

MG  =  G55  +  B25  G22  B25  ^ 

Equation  (15)  gives  the  vector  of  resistor  tree  branch  voltages  in  terms  of 
quantities  that  are  all  known.  Then,  Vj ,  I2,  and  Ig  can  be  solved  by  Equa¬ 
tion  (2),  yielding 

h  -  °22  V2  <17> 

and  I.  =  Vr  (18) 

D  *j\*  O 

respectively.  The  two  approaches  differ  in  the  size  of  the  matrix  to  be  in¬ 
verted.  The  first  approach  requires  the  inversion  of  a  matrix  (Mr)  contain¬ 
ing  the  number  of  rows  and  columns  equal  to  the  number  of  class -2  elements 
in  the  network.  The  second  approach  requires  the  inversion  of  a  matrix  (Mq) 
containing  the  number  of  rows  and  columns  equal  to  the  number  of  class -5 
elements  in  the  network.  Networks  containing  resistors  that  are  all  constant 
require  only  one  matrix  inversion.  There  is  no  practical  difference  between 
the  two  approaches.  If,  however,  the  network  contains  at  least  one  variable 
resistor,  a  matrix  must  be  inverted  at  each  solution  time  increment.  Hun¬ 
dreds  or  even  thousands  of  matrix  inversions  are  necessary,  and  the  size  of 
the  matrix  becomes  of  significant  importance.  SCEPTRE  will  automatically 
determine  which  of  the  two  approaches  should  be  taken  for  each  individual 
network. 


2.  3.  2  SOLUTION  OF  CAPACITOR  QUANTITIES 


The  capacitor  quantities  that  must  be  solved  at  each  time  step  are  Ij, 

I4,  Vj  and  V4 .  Vector  itself  will  have  been  updated  by  the  integration 
routine  and  hence  will  be  known.  The  fundamental  relationships  for  capacitors 
permit 


■  S44 


(19) 


■  S11  !1 


(20) 


Equations  (5),  (9),  (19),  and  (20)  may  be  combined  to  obtain 

ri  ’  Ms"‘  {  -B14  S44  [  B24T  '2  -  B34T  h  +  B84T  J8  ]  '  B17  B7}<21> 


where 


+  B14  S44 


(22)  j 


At  this  point,  vector  Ij  has  been  isolated  in  terms  of  all  known  quantities. 
There  remains  to  obtain  I4,  Vj  ,  and  V4  from  equations  (5),  (l),  and  (19), 
respectively. 


2.  3.  3  SOLUTION  OF  INDUCTIVE  QUANTITIES 


• 

The  inductive  quantities  that  must  be  solved  at  each  time  step  are  I3, 
V3  ,  V6,  and  Ig.  Vector  I3  will  have  been  updated  at  the  start  of  each  time 
step  and  will  be  known.  The  fundamental  relations  for  inductors  permit 


33  *3 
63  J3 


+  L36  *6 
+  L66  !6 


(23) 

(24) 


Equations  (3),  (10),  (23),  and  (24)  may  be  combined  to  obtain 
*3  =  ML  1  {  ”B34  V4  "  B35  V5  "  B37  E7  '  [B36  L66  B86 
+  L36  B86T] 


(25) 
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where 


ML  L33  f  L36  B36  +  B36  L63  +  B36  L66  B36  <26) 

Now  I3  is  written  in  terms  of  known  quantities.  Following  this,  V3 ,  V(j ,  and 
I6  may  be  obtained  from  Equations  (10)  and  (23),  (10)  and  (24),  and  (7),  respectively. 

2.  3. 4  SOLUTION  OF  VOLTAGE  SOURCE  CURRENTS  AND  CURRENT 
SOURCE  VOLTAGES 

A  complete  list  of  possible  outputs  of  a  network  would  include  the 
current  through  voltage  sources  and  the  voltage  across  current  sources. 

These  can  be  obtained  directly  from  equations  (8)  and  (4),  respectively,  since 
the  right  sides  of  these  equations  are  known  at  this  stage  of  the  computational 
sequence.  These  steps  complete  the  formal  derivation  of  all  network  currents 
and  voltages. 

2.4  SCANNING  PROCEDURE 

2.  4. 1  IDEAL  OPERATION  (SUBMATRICES  B14,  B25,  AND  B36  =  0) 

The  series  of  matrix  operations  described  in  section  2.  3  could  very 
well  be  programmed  as  they  are  to  produce  the  solution  of  the  general 
transient  analysis  problem.  All  of  the  matrix  multiplication,  addition,  etc. 
could  be  performed  at  each  time  step  to  generate  the  necessary  currents  and 
voltages.  Research  completed  during  the  study  phase  of  this  contract  showed, 
however,  that  a  very  significant  improvement  in  computer  running  time  could 
be  achieved  by  a  more  efficient  utilization  of  the  information  contained  in  the 
B  matrix. 

The  study  report*  indicated  that  various  network  voltages  and  currents 
could  be  read  or  "scanned"  directly  from  the  rows  and  columns  of  the  B 
matrix  respectively.  This  can  be  put  more  precisely  by 

VL  =  -B  VTB  (27) 

and  ITB  =  BT  II  (28) 

where  I^g,  II,  V>tb  and  Vl  represent  the  vectors  of  tree  branch  currents, 
link  currents,  tree  branch  voltages,  and  link  voltages,  respectively.  If  the 
vectors  and  the  B  matrix  are  partitioned  according  to  the  form  of  figure  2. 
equations  (27)  and  (28)  lead  to  the  first  six  equations  of  section  2.3.  Quan¬ 
tities  V2  and  I5  are  explicitly  written  in  terms  of  known  quantities  of  equa¬ 
tions  (2)  and  (6)  if  submatrix  B25  =  0.  Once  the  resistive  quantities  are 

*  AFWL-TR-65- 101,  Volume  I 
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determined,  equation  (5)  presents  I4  in  terms  of  known  quantities  if  submatrix 
B14  =  0.  In  the  same  manner  equation  (3)  presents  V3  in  terms  of  known 
quantities  if  submatrix  Bgg  =  0.  Clearly,  then,  the  condition  "B14,  B25, 
b36  =  0”  permits  solution  of  V«,  I5,  I4.  V3,  Ig,  and  Vj  directly  by  scanning 
ana  without  any  matrix  manipulation.  Quantities^,  Vg,  Ij,  Vg,  V4,  and  I3 
can  then  he  determined  from  the  operations  given  or  implied  in  the  non  zero 
portions  of  equations  (11),  (12),  (21),  (24),  (19),  and  (25),  respectively. 

Once  all  of  these  quantities  are  explicitly  obtained  in  terms  of  known  quantities, 
the  equations  are  stored,  compiled,  and  executed  at  each  time  increment 
without  recourse  to  repeated  matrix  manipulation. 


2.  4.  2  OPERATION  WHEN  B25  *  0 

In  most  large  networks,  submatrix  B25  does  not  equal  zero.  When 
this  happens,  quantities  V2  and  Ig  cannot  be  ’’scanned”  out  and  either  Equa¬ 
tion  (13)  or  (15)  must  be  solved.  Both  of  these  equations  require  the  inver¬ 
sion  of  a  matrix;  thus,  some  matrix  manipulation  must  be  done.  To  illustrate 
the  procedure,  assume  that  some  hypothetical  network  has  given  rise  to  the 
B  matrix  shown  in  figure  3.  The  network  is  fairly  typical  in  that  B14,  B3g 
-  0,  but  B25  /  0.  The  nature  of  submatrix  B25  (outlined  in  figure  3)  is  such 
that  resistors  R3  and  Rj)  contribute  no  rows  or  columns,  and  the  scanning 
process  immediately  yields 


VR3  =  E1  -  VC2  and  rRD  =  rLl>  fr0m  Which  folloW 
!R3  =  <E1  '  VC2*/,R3  and  VRD  =  ”d  :L1 


These  quantities  will  be  updated  at  each  time  step  without  matrix  manipula¬ 
tion.  The  rest  of  the  resistive  quantities  in  the  network  may  be  solved  by 
the  matrix  manipulation  implied  in  equations  (13)  and  (14).  Once  the  resis¬ 
tive  quantities  are  determined,  the  remaining  network  quantities  may  be 
scanned  out  and  stored.  The  only  repeated  matrix  manipulation  required  will 
be  for  the  three  resistor  link  currents  (Ir^,  Ir2>  IR4)  that  cannot  be  scanned 
because  Bgg  ^  0. 
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Figure  3.  B  Matrix  from  Hypothetical  Network 

2.  5  SEMI-AUTOMATIC  SOURCE  DERIVATIVES 

The  need  for  source  time  derivatives  is  established  in  equations  (9)  and 
(10)  where  E7  and  jg  are  required.  In  situations  where  non -zero  source  time 
derivatives  are  needed,  the  user  must  supply  them.  These  situations  are 
subsequently  described. 

2.  5.  1  VOLTAGE  SOURCE  IS  VARIABLE  AND  BJ7  ^  0 

If  B17  £  0  and  the  voltage  source  is  constant,  SCEPTRE  will  automat¬ 
ically  supply  a  zero  derivative.  In  the  case  where  a  non-zero  derivative  is 
required  and  the  user  fails  to  supply  it,  the  run  will  be  terminal od  with  a 
diagnostic  message.  The  situation  is  best  recognized  by  the  presence  of 
any  circuit  loop  composed  solely  of  capacitors  and  at  least  one  variable 
voltage  source. 

2.  5.  2  CURRENT  SOURCE  IS  VARIABLE  AND  B8g  ^  0 

If  Bgg  ^  0  and  the  current  source  is  constant,  SCEPTRE  will  automat¬ 
ically  supply  a  zero  derivative.  In  the  case  where  a  non-zero  derivative  is 
required  and  the  user  fails  to  supply  it,  the  run  will  be  terminated  with  a 
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diagnostic  message.  The  situation  is  best  recognized  by  the  presence  of  any 
circuit  cut  set  composed  solely  of  inductors  and  at  least  one  variable  current 
source. 

2.6  LINEARLY  DEPENDENT  SOURCES 

The  ability  of  SCEPTRE  to  correctly  process  linearly  dependent  sources 
permits  the  user  to  define  current  and  voltage  sources  that  are  linear  functions 
of  resistor  currents  and  voltages.  This  feature  is  expected  to  be  most  useful 
for,  but  not  limited  to,  applications  involving  families  of  small-signal  tran¬ 
sistor  equivalent  circuits  such  as  shown  in  figure  4. 

These  linearly  dependent  sources  require  special  treatment  because 
their  magnitudes  are  direct  functions  of  quantities  that  are  not  state  variables. 
Unless  these  sources  are  specifically  processed,  they  will  be  updated  at  the 
nth  time  step  according  to  the  values  of  their  independent  variables  at  the 
(n-1)  time  step  (as  they  would  bo  in  PREDICT,  for  example).  This  results 
in  a  "computational  delay, "  which  can  lead  to  large  errors  throughout  the 
entire  network. 


Linearly  dependent  sources  are  provided  for  in  the  general  B  matrix 
by  the  Y  andX  classification.  When  these  sources  exist,  SCEPTRE  will  set 
up  the  B  matrix  as  shown  in  figure  5.  Under  these  circumstances,  equations 
(2)  and  (6)  c  m  be  extended  to: 


V2 

'  -B24  V4  -  B25  V5  ‘  B27  E7  '  B2y  EY 

(29) 

>5 

=  B25T  52  +  B35T  ‘3  *  B85T  J8  *  Bx5T  ^ 

(30) 

Figure  4.  Low- Frequency  h -Parameter  Equivalent  Circuit 
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Figure  5.  B  Matrix  with  Linearly  Dependent  Sources 

Since  any  resistor-voltage -dependent  voltage  source  must  depend  on  a  resis¬ 
tor  tree  branch  voltage  or  a  resistor  link  voltage,  there  arises 

EY  =  kj  V2  +  k2  V5  (31) 

where; 

•  kj  is  a  matrix  of  constants  containing  the  number  of  rows 
equal  to  the  number  of  class -Y  voltage  sources  and  the 
number  of  columns  equal  to  the  number  of  nonscannable 
class-2  elements. 

•  kg  is  a  matrix  of  constants  containing  the  number  of  rows 
equal  to  the  number  of  class -Y  voltage  sources  and  the 
number  of  columns  equal  to  the  number  of  nonscannable 
class-5  elements. 

•  EY  is  a  vector  composed  of  class -Y  voltage  sources. 

Also,  since  any  resistor-current -dependent  current  source  must  depend  on 
a  resistor  tree  branch  current  or  a  resistor  link  current,  there  arises 

JX  =  k3  I2  +  k4  Ig  (32) 

where: 

•  kg  is  a  matrix  of  constants  containing  the  number  of  rows 
equal  to  the  number  of  class-X  current  sources  and  the 
number  of  columns  equal  to  the  number  of  nonscannable 
class-2  elements 
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•  k4  is  a  matrix  of  constants  containing  the  number  of  rows 
equal  to  the  number  of  class -X  current  sources  and  the  number 
of  columns  equal  to  the  number  of  nonscannable  class -5 
elements 

•  JX  is  a  vector  composed  of  class-X  current  sources. 

If  Equations  (31)  and  (32).  together  with  V2  =  R22I2  and  *5  "  ^5  Mi 
are  substituted  into  Equations  (29)  and  (30) 

R22  l2  +  D25  V5  +  B2y  [kl  R22  *2  +  k2  V5  ]  '  "B24  V4  ’  B27  E7  (33^ 

G55  V5  "  B25  l2  '  Bx5  [k3  X2  "  k4  G55  V5]  “  B35  X3  +  B85  J8 


The  last  two  equations  may  be  consolidated  as 


(R22  +  B2y  kl  R22^  ^B25  +  B2y  k2^ 

V 

~B24  V4  '  B27  E7 

.‘-B25T-Bx5Tk3><G55-Bx5Tk4G55>  . 

B35  X3  +  B85  J8 

If  the  large  matrix  on  the  left  side  is  called  MRG,  then 


-1 

”B24  V4  ‘  B27  E7 

=  MRG 

„  T  .  _  T 

B35  X3  +  B85  J8 

and  the  resistive  quantities  I2  and  V5  can  be  determined  without  computation¬ 
al  delay.  Note  that  since  all  four  "k”  matrices  are  constrained  to  be  constant, 
the  linearly  dependent  source  feature  itself  will  not  require  any  more  than 
one  inversion  of  MRG.  If  any  variable  resistors  are  present  in  a  network, 
MRG  must  of  course  be  inverted  at  each  time  step.  In  addition,  the  extra 
row  and  column  that  was  added  to  the  D  matrix  of  figure  2  by  the  linearly 
dependent  sources  will  add  terms  to  the  equations  used  in  the  solution  of 
capacitive,  inductive  and  some  source  quantities.  Specifically,  equations 
(5),  (3).  (4)  and  (8)  become 

'4  =  BHT  h  +  B24T  !2  +  B34T  *3  +  B84T  J8  +  Bx4T  ^  <5'> 
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V3  a  “B34  V4  -  B35  V5  *  B36  V6  *  B37  E7  '  B3y  W 


8 


B84  V4  -  B85  V5  *  B86  V6  '  B87  E7  “  %  EY 


I»  =  B,  „T  I.  +  BnJ  U  +  Ba„T  Ig  +  B  aJ  J0  4  Bl  JX 


‘7  "  17  ‘1  T  27  *2 T  “37 

and  the  new  relations 


'87  ”8  x7 


m 

(8*) 


™  *  -°*4  V4  -  Bx5  V5  -  Bx7  E7  *  Bxy  EY 

(35) 

'Y=VVB3yVB84TVBx> 

(36) 

now  exist. 

A  restriction  must  be  placed  on  these  sources  based  on  the  content  of 
Section  2.  5.  The  B  matrix  o!  figure  5  transforms  equations  (9)  and  (10) 
into: 


V1  =  -B14V4  Bl7E7-BlyEY 
j6  =B36Ti3*B86TVB*9TjX 


00*) 


Now,  additional  time  derivatives  EY  and  JX  are  required  whenever  B*„  or 
Bx6  *  »•  ’ 


Differentiation  of  equations  (31)  and  (32)  would  involve  quantities  V2 .  V5 
i2,  and  {5.  The  formulation  contains  no  provisions  for  these  quantities  and  * 
there  is  no  way  the  user  could  know  them  to  supply  them  as  input  data. 
SCEPTRE  will  automatically  check  for  the  existence  of  non- aero  Biyor  Bx6 
and  terminate  the  run  with  a  diagnostic  message  when  they  occur.  The  situ¬ 
ation  can  be  recognized  by  the  presence  of  any  circuit  loop  composed  solely 
of  capacitors  and  at  least  one  linearly  dependent  voltage  source,  or  the  pres¬ 
ence  of  any  circuit  cut  set  composed  solely  of  inductors  and  at  least  one 
linearly  dependent  current  source. 

2.  7  INITIAL  CONDITIONS  SOLUTION 

2.  7.  1  TECHNIQUE  DESCRIPTION 

Many  practical  circuits  require  computer  solution  of  the  Initial  condi¬ 
tions  prevailing  at  the  start  of  the  transient  (time  -  lQ)  before  the  transient 
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solution  can  begin.  These  values  can  always  be  determined  by  a  separate 
transient  run  in  which  all  forcing  functions  are  held  at  the  values  for  time 
-  tQ.  However,  an  alternate  procedure  based  on  an  iteration  technique  using 
independent  variables  other  than  time  was  considered  desirable.  This  pro¬ 
cedure  presents  the  advantage  of  economy  of  machine  time  on  all  circuits  for 
which  convergence  occurs.  This  section  will  describe  the  formulation  of  this 
portion  of  SCEPTRE,  which  is  completely  independent  of  the  transient 
formulation. 

If  an  L  tree  is  set  up  based  on  the  preference  order  E,  L.  R,  C,  a 
general  B  matrix  may  be  set  up  according  to  the  procedure  outlined  in 
section  2.  2  The  zero-valued  submatrices  arise  from  the  L  tree  and  the 
preference  order*.  The  resulting  B  matrix  is  shown  in  table  n. 

Table  II 

GENERAL  B  MATRIX  FROM  THE  L  TREE 


4 

5 

6 

1 

7  1 

1 

B14 

B15 

B16 

Bn 

2 

0 

B25 

B26 

B27 

3 

0 

0 

B36 

B37 

00 

B84 

B85 

B86 

B87 

05 

B94 

B95 

B96 

B97 

0 

B04 

B05 

B06 

B07 

The  following  equations  (among  others)  arise  from  this  B  matrix  if 
vectors  V$,  I4  and  submatrix  B94  are  assumed  to  be  zero.  These  assump¬ 
tions  are  based  on  the  known  final  values  of  and  I4  for  the  initial  condition 


♦For  example,  the  submatrix  B24  must  always  be  zero  since  non -zero  entries 
in  it  could  only  arise  when  resistor  iinks  close  loops  containing  capacitor  tree 
branches.  The  preference  order  prohibits  this  possibility. 
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problem  and  the  absence  of  any  current-source  capacitor  cut  sets  (see  the 
restrictions  in  section  2.7.4). 


V9  +  B95  V5  +  B97  E7  =  0 

(37) 

V2  4  °25  V5  +  D27  E7  ~  0 

(38) 

X5  "  D25  l2  "  B85  J8  “  B95  J9  “  B05  J0  =  ° 

(39') 

To  get  all  the  variables  of  equation  (39')  in  terms  of  Vg,  V2 ,  V5  ,  and 
effective  sources,  the  following  substitutions  are  made: 


C55  V5 


G22  V2 


Jo  "  Qj< 


J9  G99  V9  +  ^9 


where: 

•  a  is  a  matrix  containing  the  number  of  rows  equal  to  the 
number  of  secondary  current  sources  and  the  number  of 
columns  equal  to  the  number  of  primary  current  sources 

•  G55  and  G22  are  diagonal  matrices  containing  only 
conductances 

•  Ggg  is  a  diagonal  matrix  containing  only  diode  and  transistor 
junction  conductances 

•  Q  terms  are  described  in  appendix  II. 


Then,  equation  (39')  becomes 


G55  V5  "  B25  G22  V2  "  D85  J8 


1 

1 

D95  +  °05  a 

°99  V9  '  % 

—  — 

(39) 
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where  the  Jacobian  is 


If  . 


dF 

dv‘  <V  V  V 

dF 

-W~  <v9'  V2‘  V5> 

dF, 

Ws  (V9-  V2-  V9* 

avj  (V9’  V2'  V5> 

dF,  . 

dV2  <  V9-  V2'  V 

dr0 

v2-  V 

dF« 

av9  < v9'  V2.  V5> 

dF, 

av2  (v9-  v2-  v5> 

dF, 

-d-vf  (V9'  V2-  V 

(41) 

AVg  -  Vg  (n  -  1)  -  Vg  (n) 

AV2  V2  (n  +  X)  "  V2  (n) 

AVj-  V5  (n  +  1)  -  V5  (n) 

where  n  is  used  to  designate  the  results  of  the  nlh  iteration  pass,  then 
equation  (40)  becomes 


F1  (V9’  V2’  V5) 

Vg  (n  4  1)  -  Vg  (n) 

F2  <V9’  V2’  V5> 

+  z 

V0  (n  +  1)  -  V2  (n) 

F3  <V9'  V2’  V5> 

V5  (n  •  D  -  V5  (n) 

_  - 

leading  directly  to 


Vg  (n  *  1) 

Vg  (n) 

F,  (V9.  V,.  V5> 

v2  0.  1) 

v2  (,.) 

-  y._1 

F2  (Vy.  V2.  Vy) 

jv5  (n  ■  1) 

V,  (n) 

L  0 

'  (V9-  V2-  V 

* 

(42) 


(4:t) 
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Equation  (41 )  may  be  written  more  explicitly  if  the  indicated  differentiations 
are  performed  on  equations  (37),  (38),  and  (39)  to  obtain: 


Z  = 


O 

I 


B 


95 


B 
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-  [B95T  *  B05T  °] 
so  that  equation  (43)  becomes 


B25  G22  °55 


-1 


Vg(n "  1) 

V"> 

-  1 

Fi 

(Vg,  V2'  V5> 

V2(n  +  1) 

= 

v2<"> 

- 

z 

F2 

<V9-V2-V5> 

V5(n  +  1> 

V5<n> 

F3 

(Vg-  V2>  V 

L.  - 

Equation  (44)  may  be  written  in  more  convenient  form  as  follows  (see  also 
Appendix  IV): 


(44) 


vg  (n  -  1) 

— 

I 

O 

B95 

V2  (n  ♦  1) 

O 

I 

B25 

V5  (n  .  1) 

-[B9ST'B0ST°]S9 

"B25  G22 

C55 

B97  E7 


?95T  +  B05T  a]  Q9  +  B85  J8 


Equation  (44*)  signifies  a  computational  sequence  as  follows: 
Quantities  Vg,  V2,  and  V5,  the  vectors  of  voltages  across  primary  current 
sources,  resistor  link,  and  resistor  tree  branches  respectively,  each  have 
some  assumed  value  (usually  zero)  to  begin  the  computation  at  n  =  0.  AU 
members  of  the  right  side  of  equation  (44')  that  may  be  voltage  dependent 
are  updated.  The  left  side  of  equation  (44')  is  then  computed  and  the  first 
iteration  (n=l)  is  complete.  The  right  side  of  equation  (44*)  Is  re-evaluated 
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on  the  basis  of  the  results  of  the  first  iteration  and  the  left  side  is  again  com¬ 
puted,  thereby  completing  the  second  iteration  (n=2).  This  process  is  repeated 
up  to  100  times.  After  any  Iteration,  if 

lV(„.  1)-V(n,  I*0-' 001  I  V.  1)1  «  0001  <45> 

is  satisfied  for  the  set  of  all  voltages  in  equation  (44),  convergence  is  con¬ 
sidered  to  have  occurred  and  the  procedure  is  terminated.  If  after  100 
iterations  equation  (45)  is  not  satisfied,  a  diagnostic  message  will  be  printed 
indicating  that  convergence  has  not  occurred.  Experience  to  date  has 
shown  that  the  Newton-Raphson  procedure  will  converge  for  most  circuits  in 
loss  that  30  iterations. 

Convergence  of  the  Newton-Raphson  method  can  sometimes  be  pro¬ 
hibitively  delayed  if  for  some  reason  a  large  forward  bias  (V>  0. 8v)  is  applied 
to  any  diode  or  transistor  junction  that  has  been  represented  by  the  user  in 
the  conventional  closed  form  J  =  Ig(e®V_i),  jn  this  case,  the  slope  of  the 

rtliMip  curve,  given  by  Oj  »  *Jy  will  ronlrlbutr  a  very  large  term 

in  Ggg  and  consequently  in  the  Z  matrix.  The  practical  results  of  this  can  be 
more  easily  appreciated  by  consideration  of  a  one  equation  system  as  repre¬ 
sented  by  the  second  equation  in  Appendix  HI.  A  large  derivative  caused  by  a 
highly  forward-biased  diode  leads  to  a  very  small  step  (AX)  in  the  indepen¬ 
dent  variable.  Many  steps  will  therefore  be  required  to  complete  convergence. 
This  situation  is  avoided  in  SCEPTRE  by  the  inclusion  of  a  subroutine  called 
DISCLF.  This  subroutine  effectively  ensures  that  the  true  operating  point  on 
the  diode  curve  is  approached  in  the  iteration  process  from  a  lower  voltage 
rather  than  a  higher  voltage.  Smaller  diode  slopes  are  used  so  that  larger 
(and  fewer)  steps  can  be  taken  toward  convergence.  This  is  described  in 
reference  2. 


If  convergence  has  occurred,  quantities  Vg,  Vg,  and  V5  arc  now  known. 
There  remains  to  compute  only  capacitor  voltages  and  inductor  currents  since 
capacitor  currents  and  inductor  voltages  must  be  zero. 

2.  7.2  COMPUTING  CAPACITOR  VOLTAGES 

From  the  B  matrix  in  table  n  the  capacitor  link  voltages  can  be  written 
in  terms  of  the  tree  branch  voltages  as 


-BH  V4  -  °15  V5  -  D17  E7 


(4I>) 


In  addition  the  principle  of  conservation  of  charge  permits 

-11,4  r  Cn  V,  .  C44  V4  -Di4T  Cj ,  V,  (0)  •  C44  V4  (0)  (47) 

where  Vj  (0)  and  (0)  are  initial  voltages  that  may  be  specified  by  the  user. 
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Equations  (46)  and  (47)  lead  to 

V4  -  Mc  '*  {-BHT  C11  B1S  V5  -  B14T  C11  Bn  Sj  '  B14T  CU  V1  (°> 

+  C44  V4  (0) }  <4B> 

T 

where  M  =  B.,  Cn  B..  +  CAA 
c  14  11  14  44 

Cnee  V4  is  determined,  Vi  may  be  found  from  equation  (46). 

Note  that  Class-4  elements  can  occur  only  if  a  capacitor  cut  set  exists.  ^ 
Otherwise  equation  (48)  will  be  identically  zero  and  all  capacitor  voltages  can 
be  determined  from  equation  (46). 

The  use  of  equation  (47)  permits  solution  of  a  class  of  networks  in 
which  the  final  value  of  capacitor  voltage  is  dependent  on  an  initial  value. 
Consider  the  circuit  shown  in  figure  6.  which  contains  a  capacitor  cut  set. 


Figure  6.  Capacitor  Cut  Set  Circuit 

If  both  capacitors  are  initially  uncharged,  the  final  values  of  the  cap  .tor 
voltages  after  the  switch  is  closed  must  be  VC1  -  1  volt,  VC2  -  9  volts.  If 
however,  capacitor  Cl  has  an  initial  charge  of  5  volts,  the  principle  of  con¬ 
servation  of  charge  (reflected  in  equation  (47))  requires  the  final  result  to 
be  VC1  =  5.  5  volts,  VC2  =  4.  5  volts. 

2.  7.  3  COMPUTING  INDUCTOR  CURRENTS 

Since  V2  is  known, 

>2  '  B22_1  V2  <49> 

The  inductor  tree  branch  currents  then  can  be  written  from  the  B  matrix 
in  terms  of  link  currents  as 


*6  B26 


T  T 

l2  +  B36  k 


B 


86 


V  B96 


T  T 

J9+B06  JC 


(50) 
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In  addition,  the  flux  relations  around  an  inductor  loop  permit 

B36  L66  *6  +  L33  *3  =  B36  L66  *6  (0)  +  L33  *3  (0)  (5l) 

where  Ig  (0)  and  I3  (0)  are  initial  currents  that  may  be  specified  by  the  user. 
Equations  (50)  and  (51)  lead  to 


-1 


*3  =  ML  {  L33  *3  (0)  +  B36  L66  *6  (0) 


’b36  l66 


Brt„T  I.  +  B0„TJ0"  J9  +  Bq6T  JQ1  \ 


26  2  86  8  96 


(52) 


where 


ML  "  B36  L66  B36  +  L33 


Once  I3  is  determined.  Ig  may  be  found  from  equation  (50).  Note  that  Class-3 
elements  can  occur  only  if  an  inductor  loop  exists.*  Otherwise,  equation 
(52)  will  be  identically  zero  and  all  inductor  currents  can  be  determined  from 
equation  (50). 

The  use  of  equation  (52)  permits  the  solution  of  a  class  of  networks  in 
which  the  final  values  of  inductor  currents  are  dependent  on  the  sizes  of  the 
respective  inductances.  Consider  the  circuit  shown  in  figure  7,  which  con¬ 
tains  an  inductor  loop.  If  both  inductors  are  initially  relaxed,  the  final  values 
of  the  inductor  currents  after  the  switch  is  closed  must  be  IL1  =  1  amp. 

IL2  =  9  amps. 


L2  -  lh 


Figure  7.  Inductor  Loop  Circuit 


"See  paragraph  (subsequent  section  2.7.4)  on  network  restrictions  for  quali¬ 
fication  of  this. 
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The  task  of  computing  the  capacitor  voltages  and  inductor  currents  com¬ 
plete  the  initial -conditions  problem.  These  quantities  may  then  be  trans¬ 
ferred  to  the  transient  program  to  serve  as  initial  conditions. 

2.7.4  RESTRICTIONS 

Since  the  initial-conditions  program  is  formulated  differently  and  for  a 
different  purpose  than  is  the  transient  program,  certain  restrictions  apply 
to  the  initial-conditions  program  that  do  not  apply  to  the  transient  program. 
These  restrictions  and  the  practical  considerations  that  lead  to  them  are  as 
follows: 

1.  No  circuit  containing  a  loop  composed  entirely  of  voltage 
sources  and  inductors  will  be  accommodated.  This  situation 
would  cause  an  infinite  inductor  current  and  is  obviously  of 
no  practical  importance.  The  presence  of  an  E-L  loop  is 
disclosed  by  the  condition  B^  *  0. 

2.  No  circuit  containing  a  cut  set  composed  entirely  of  current 
sources  and  capacitors  will  be  accommodated  (see  figure  8). 

This  situation  would  invalidate  equation  (37)  and  complicate 
the  solution  process  by  requiring  that  V4  be  carried  along  in 
the  Newton- Rap hs on  procedure.  This  cut  set  situation  can 
always  be  removed  by  arbitrarily  connecting  a  large  resistor 
from  note  A  to  the  ground.  Note  that  the  configuration  of 
figure  8  could  be  handled  by  the  transient  portion  of  SCEPTRE 
if  a  Newton-Raphson  solution  was  not  desired.  The  pres¬ 
ence  of  a  J— C  cut  set  is  disclosed  by  the  condition  that 
either  B  34.  B94.  or  B04  t  0. 


- ^ - 
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>  —  —  ' 

Figure  8.  A  Current -Source  Capacitor  Cut  Set 
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3.  The  choice  of  independent  va r  tables  will  bo  .somowluii  r*  si  ni  :«•<!. 

No  resistor  or  inductor  current  may  bo  usod  as  an  independent 
variable.  Furthormoro,  a  capacitor  voltage  may  bo  usod  as  an 
indopondont  variable  only  if  it  is  in  parallel  with  a  resistor  or 
current  source  (Vc  =  Vr,  Vc  Vj).  The  objective  hero  is  to  avoid 
the  need  for  any  auxiliary  computation  between  passes  of  the 
iteration  procedure. 

4.  Networks  containing  entirely  capacitor  cut  sets  can  bo  accommo¬ 
dated  only  if  the  members  of  the  cut  set  are  constant. 

5.  Networks  containing  only  inductor  loops  can  be  accommodated  only 
if  the  members  of  the  loop  are  constant. 

2.  8  INTEGRATION  ROUTINE 


Three  integration  routines  will  be  optionally  available  for  use  in 
SCEPTRE.  Two  of  these.  RUK  and  TRAP,  were  available  in  PREDICT 
and  have  been  only  slightly  modified.  The  third  routine,  called  XPO,  was 
developed  at  the  IBM  Scientific  Center  in  Palo  Alto,  California,  by  Dr.  R. 
Wartcn  and  Mr.  M.  Fowler  and  adapted  for  use  in  SCEPTRE.  Studies 
to  date  indicate  that  XPO  is  usually,  although  not  always,  faster  than  the 
other  methods.  For  that  reason,  this  routine  will  always  be  used  unless 
the  user  explicitly  requests  otherwise  in  the  RUN  CONTROL  section  of  any- 
run. 

2.  8.  1  RUK  FORMULAS  AND  VARIABLE  STEP  SIZE  CONTROL 

The  well-known  Runge-Kutta  fourth-order-accurucy  formulas  (Ref.  2) 
for  the  numerical  solution  of  the  differential  equation 

y  f  (t.  y) 

are  uivon  by 

v  (i  •  h)  y  (l)  •  1/0 [kj  -  2k2  *  2k;j  .  k4  ]  (r>3) 

where 

k,  h  f  [t,  y  (I)] 

k2  '•  1  |  1  >  2  *  v  0)  ’  2kl| 
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Now  let 

ur  u2’  *  r  *2  ~  0 

Set 

Uk  ‘  U1  +  u2  I  yk  (t  +  h) 

Lk  =  *1  T  *2  I  yk  +  h) 

If  |  |  >  1.  5  Ufc  for  some  k,  the  integration  step  h  is  halved,  the  in¬ 

dependent  variable  is  restored  from  t  +  h  to  t,  and  the  values  of  y  and  y  are 
restored  to  the  values  at  time  t. 

If  |  e^  |  >  0.  75  for  some  k,  and  |  ejJ  <  1.  5  Uu  for  all  k,  the  current 
integration  step  is  accepted,  but  the  step  size  is  halved  for  succeeding  steps. 

If  <  I  e^  I  <  0.  75  Ur  for  all  k,  the  step  size  is  unaltered. 
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If  |ekl  <  Uk  for  k,  and  |ek|  <  Lkfor  at  least  one  k,  a  doubling  in¬ 
dicator  is  activated.  Actual  doubling  is  delayed  for  seven  time  steps.  Halving 
always  takes  precedence  over  doubling;  thus  anytime  a  halving  signal  is  re¬ 
ceived,  the  step  size  is  halved  and  doubling  is  delayed  for  at  least  seven  steps. 
Similarly,  after  successful  doubling,  another  seven  steps  must  elapse  before 
the  step  size  can  be  doubled  again. 

Recommended  choices  for  uj,  U2,  ll ,  and  I2  are  as  follows.  If 
absolute  error  control  is  desired, 

set 


Uj  =  0.0075  u2  =  0 

1  =  0.00005  12  =  0. 

If  relative  error  control  is  desired, 
set 

Uj  =  0.005  u2  =  0.005 

lx  =  0.00005  12  =  0.00005. 

If  smaller  step  sizes  are  desired  than  the  ones  yielded  by  the  above 
settings,  the  0.  0075  and  0. 00005  settings  should  be  reduced  by  a  factor  of  32. 
This  will  yield  half  the  previous  step  sizes. 

2.  8.  2  MODIFIED  TRAPEZOIDAL  INTEGRATION  (TRAP) 

2.  8.  2. 1  Method 

Given  the  differential  equation 

y  (t)  =  f  (t,  y  (t»,  y  (0)  =  yQ. 

Consider  the  following  numerical  integration  scheme, 

yp  (l  +  !  h)  =  y  (t)  +  \  y  (I)  <r,r,> 
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yp  (t  +  |h)  =  f[t  +  |h,  yp(t+|h)] 

-  f  [t  +  |  h,  y  (t)  +  £  y  (t)] 

(56) 

yc  (t  +  h)  =  y  (t)  +  J  [2  y  (t)  +  yp  (t  +  |  h)] 

(57) 

yc  (t  +  h)  =  f  [t  +  h,  yc  (t  4  h)j 

(58) 

Step -size  control  is  achieved  by  computing 


E  =  §h 


vp  (t  +  |  h)  -  y  (t) 


(59) 


The  magnitude  of 


indicates  any  changes  to  be  made  in  the  step  size. 

2.  8.  2.  2  Truncation  Error 

Inasmuch  as  y  =  +  y  fy  and  the  true  solution  y^,  is  approximated  by 

.2 

yT  (t  +  h)  =  y  (t)  +  h  y  (t)  +  j  y  (t). 


one  obtains, 

.2 

yT  (t  +  h)  =  y  (t)  +  h  y  (t)  +  -j-  (ft  *  y  fy). 

On  the  other  hand 

f  [t  +  y  ,  y  (t)  +  |  y  (t)]  =  f  (t,  y  (t)>  +  y  ft  +  y  fy 
whence  Equation  (57)  becomes 

y c  (t  +  h)  =  y  (t)  +  j  |^2y  (t)  +  y  (t)  +  -g-  f^  +  ^  Y  fy  j 
Subtracting  Equation  (57)  from  (60)  yields 


yT  (t  +  h)  -  yc  (t  +  h)  = 


(60) 


(61) 


30 


which  in  turn  yields  the  approximate  local  truncation  error. 

From  equations  (56)  and  (61)  one  gets 

•  /.  3h\  */x\*3h.  h  •  , 
yp  (t  +  y )  -  y  (t)  =  y  ft  +  2  y  f 

whence 

2 

2h  r  •  ,,  3h  ,  •  /.,1  .  ,2.  h  •  , 

y  [  y  (t  +  y  )  -  y  (t)  J  =  h  f(  +  j  y  iy  (62) 

One  notes  that  if  ft  =  0,  then  equation  (62)  is  a  good  representation  of 
the  local  truncation  error,  neglecting  higher  order  terms. 

Experience  indicates  that  for  systems  of  equations  of  the  form  Y  = 

A  (Y)  Y  +  B,  where  A  is  piecewise  constant,  the  method  as  well  as  step-size 
control  is  adequate. 

2.  8.  2.  3  Stability 

Consider  the  differential  equation 

y  =  -ay,  y(0)  =  yQ,  a>0 


The  numerical  integration  scheme  given  by  equations  (55)  through  (58) 

yields 


This  is  easily  established  by  induction.  Thus,  for  numerical  stability 
it  is  necessary  and  sufficient  that  j  1  -  ah  +  j  <  This  yields  0  < 

ah  <  6  as  shown  in  Figure  9. 


Inasmuch  as  |  ah  |  <  6  is  necessary  and  sufficient  for  numerical  stability, 
we  say  that  the  modified  trapezoidal  method  has  a  stability  radius,  r,  equal  to 
6.  Moreover,  from  equations  (55)  and -(56)  we  see  that  the  method  requires 
two  derivative  evaluations  (passes)  per  integration  step.  Thus,  the  pass  num¬ 
ber  p  associated  with  the  method  is  2. 


The  ratio  r/p  is  a  measure  of  a  method's  efficiency  in  (ho  sense  of 
minimum  number  of  integration  steps  per  given  time  interval. 
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1  3  5  6  oh 

-1  - 

Figure  y.  Root  Locus  of  TRAP  Characteristic  Equation 
The  higher  r/p,  the  fewer  steps  are  required  for  a  given  a. 

The  following  list  shows  r,  p  and  r/p  for  a  few  representative  methods: 


METHOD 

r 

P 

r/p 

Modified  Trap. 

6 

2 

3 

Euler  (Ref.  4) 

2 

1 

2 

Trapezoidal  (Ref.  5) 

2 

2 

1 

NIDE  (1PASS)  (Ref.  6) 

0.  8 

1 

0.  8 

Runge-Kutta  (Fourth-order) 

2.  78 

4 

0.  7 

MEDE  (2  PASS) 

1 

2 

0.  5 

Hamming  (Ref.  7) 

0.  85 

2 

0.425 
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2.  8.  2. 4  Step -Size  Control 

As  mentioned  previously,  the  quantity 

is  used  for  step -size  control. 

(1)  The  method  is  essentially  the  same  as  for  RUK.  Let  Uj,  U2,  lj, 
1 2  be  real  numbers  >0. 

Compute  U  =  Uj  +  Ug  yc  (t  +  h)  j  ,  L  =  1^  +  yfi  (t  +  h)  | 

If  |  Eg  >  0.  75  U,  the  step  size  is  reduced. 

If  |  Eg  |  <  L,  the  step  size  is  increased. 

If  L  <  |  Eg  <  0.  75  U,  the  step  size  is  not  altered. 

(2)  Choice  of  Uj,  Ug,  1^,  lg. 

Recall  that  E2  attempts  to  represent  the  local  truncation  error.  For 
absolute  error  control  set 

ux  =  0.01  ,  Ug  =  0 

1  =  0. 0005  ,  lg  =  0 

For  relative  error  control  set 

Uj  =  0.  001  ,  Ug  =  0.  01 

1  -  0. 00005  ,  12  =  0. 0005 

The  above  values  work  well  in  practice. 

Reducing  the  above  settings  by  a  factor  of  4  will  reduce  the  step  size 
by  a  factor  of  2. 
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2.  8.  3  EXPONENTIAL  INTEGRATION  (XPO) 
2.  8.3.1  Method 


Given  the  differential  equation 

y  (t)  --  f  [t,  y  (t)j ,  y  (0)  =  y0, 

this  integration  routine  gives  the  computed  solution 

pXh_  * 

y  (t  +  h)  -  y  (t)  +  yA  (t)  h  +  —  h  yp  (t)  (64) 

except  for  two  cases  that  will  be  covered  later.  In  equation  (64) 

y  (t)  -  y  (t  -  hQ) 

YA  (0= - F - ~ 

no 

where  t  -hQ  is  the  last  point  computed; 


yp  (t)  =  y  (t)  -yA  (t): 

(65) 

X  =  y  (t)  /  yp  (t) 

(66) 

Since  v  (t)  cannot  be  computed  explicitly  by  the  program,  we  take  a  small 
Euler  integration  step  and  approximate  'y  (t)  by  y'g  (t): 


yg  (t  +8)  =  y  (t)  +  8y  (t) 
yg  (t  +8)  =  f  [t  +  8  ,  y&  (t  +8) 


h  <*>  = 


ys  (t  4-8)  -  y  (t) 

S 


0<S<x 


(67) 


The  first  exception  occurs  if  yg  (t)  and  yp  (t)  have  the  same  sign. 
Then  set 

YA  (0  =  °- 

which  results  in  equations  (65)  and  (66)  becoming 


Yp  (0  "  y  (t)  (68) 

x=  y8  (t)  /  y  (t). 
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This  provides  a  belter  approximation  to  the  solution  of  the  equation  y  Xy  *  b 
without  decreasing  the  overall  effectiveness  of  the  method. 

The  second  exception  occurs  if  X,  computed  either  by  equations  (66) 
or  (68)  is  non-negative.  Then  the  term 


eXh  _i 
Th 


Xh 


in  Equation  (64)  is  replaced  by  1  +  ,  thus  changing  equation  (64)  to 


y  (t  +  h)  =  v  (t)  +  yA  (t)  h  +  (l  +  tt-  )  hy  (t) 


(69) 


It  can  be  shown  with  some  effort  that  equation  (69)  is  equivalent  to  a  trape¬ 
zoidal  type  integration  method. 

2.  8.  3.  2  Truncation  Error 


Noting  that 


ex-l 


,  when  expanded  in  a  Taylor  series,  becomes 


i 

Equation  (64)  can  be  written 

y  (t  +  h)  =  y  (t)  +  yA  (t)  h  +  yp  (t)  h  +  Xyp  (t)  +  h  -  +  O  (h4) 

(70) 

=  y  (0  +  y  (t)  h  +  ^  y'g  (t)  +  ^  Xyg  (t)  +0  (h4) 


Referring  to  equation  (67),  note  that 

VS  (t)  =  y  (t)  +  y  [y  'W  -y  (0  ]  +  °  ( 32) 


(71) 
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Substituting  equatio'  (71)  into  equation  (70)  yields 


•  •  h2 

y  (t  +  h)  =  y  (t)  +  y  (t)  h  +  y  (t)  y- 


?2S  ['y  (t)  -V  (t)  ]  +  O  (h3) 


(72) 


It  is  noted  that  the  use  of  either  exception  to  the  principal  equation  (64)  still 
yields  the  truncated  expression,  equation  (72). 

Since  8  <  the  method  is  second  order  exact,  i.  e. ,  the  error  is  of 

3  ’ 

order  h  .  This  characteristic  will  be  used  in  equation  (73)  in  the  next 
section. 


2 .  8.3.3  Step- Size  Control 


For  t<t  +  £  <  t  +  h,  let  y<r  (t  •*-  £  )  be  the  true  solution  to  the  differential 
equation,  let  yc  (t  +  £  )  be  the  computed  solution,  let  ye  (t  i-  £  )  be  an  ex¬ 
pression  obtained  by  differentiation  of  yc  (t  +  £  ),  i.  e. , 


yA  (t)  +  (1  +  X£  )  yp  (t),  X  >0 
>'a  (t )  +  e  X £ yp  (t )  ,  \<  0 


and  let  yc  (t  +■  £  )  be  obtained  by  substituting  yc  (t  +  £ )  into  the  differential 
equation. 


Ideally,  step  control  should  be  based  on  the  expression 


■hr. 


Eo  =  J  q  [y  t  ^  +  £ )  -ye  +  £  )  d  c 


however,  y  (t  +£)  is  not  available  except  at  £  =  0. 


Making  use  of  the  fact  that 

yc  (t  +£  )  =  f  {t  +  £  ,  yT  (t  +  £  )  +  (yc  (t  +£ )  -yT  (t  -•-£ )]  J 
=  yx  (t  +£)  +  °  (£3) 


(73) 
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we  obtain 


Eo  -  f  [yc  (t  +  C)  -ye  (l  )]  d  £  +  /  [W  (t  ->'c  (t  +  £ )]  d  * 

=  J  [yc  (t  +  £  )  -ye  (t  ♦-  C  )]  d  c  +  0  (h4) 

If  yc  (t  +  £  )  -ye  (t  +  O  does  not  change  sign  for  0<£<h, 

|  Eq|  <h|yT(t  +  h)  -ye  (t  +  h)|<h|yc  (t  -r  h)  -ye  (t  +  h)j  +  |o  (h  )  j 

Tims  tor  small  h, 

E  -  h  |  yc  (t  ►  h)  -ye  (t  +  h)  | 
yields  an  estimate  of  the  truncation  error. 

Letui,  U2,  If,  1 2 be  real  numbers  >0.  Compute 

U  =  Uf  +  u  2 1  y  c  <t+h)  I 
L  =  li+  12  |yc  (t+h)| 

If  |  E  |  >U,  the  step  size  is  reduced.  If  |  E  j  <  L,  the  step  size  is  increased. 
If  L  <  E  <  U,  the  step  size  is  not  altered. 

For  absolute  error  control  set 
u  j  =  0.  0075  U2  -  0 

1 1  -  0.0002  1 2  =  0 

For  relative  error  control  set 

u  i  =  0.  005  u  2  =  005 

lj  -  0.0001  1 2  --  0.0002 
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2.  8.  4  GENERAL  ABSOLUTE  AND  RELATIVE  ERROR  CRITERIA 

All  three  numerical  integration  methods  described  in  this  section  are 
approximate  methods  that  are  used  to  integrate  differential  equations  and 
thereby  update  the  state  variables  at  each  time  step.  Each  of  these  methods 
has  an  automatic  control  that  allows  step  size  to  be  increased  and  decreased 
during  the  transient  problem.  The  method  of  control  used  compares  some 
function  of  the  step  size  and  the  derivatives  of  the  state  variables  to  a 
quantity  that  serves  as  a  standard  as.  for  example,  in  equation  (59).  The 
standard  may  be  a  constant  or  a  variable  function  of  the  state  variable.  If 
the  former  is  used,  it  is  termed  an  absolute  error  criterion;  if  the  latter  is 
used,  it  is  termed  a  relative  error  criterion. 

Consider  the  situation  in  which  two  state  variables  are  very  different 
in  size.  Let  yj  (t)  =  1  and  y2  (t)  =  100.  For  all  practical  purposes,  it  is 
usually  unnecessary  to  integrate  the  derivatives  of  these  state  variables  to 
the  same  accuracy.  If  an  absolute  error  criterion  is  used,  this  is  just  what 
is  done,  and  the  step  size  may  be  unnecessarily  inhibited.  If,  however,  a 
relative  error  criteria  is  used,  effectively  a  larger  error  will  be  tolerated 
for  the  integration  of  y£  (t)  and  a  larger  step  size  will  be  permitted.  In 
general,  then,  it  would  be  best  for  the  user  to  use  relative  error  control 
when  large  values  of  state  variables  (capacitor  voltages  and  inductor  cur¬ 
rents)  are  expected. 

Relative  error  control  is  programmed  with  all  three  integration  methods, 
but  the  user  may  easily  modify  the  relative  controls  or  enter  absolute  controls 
( section  2.  2.  9  of  Volume  I) .  If  larger  numbers  are  used  for  the  u^  or  U2 
entries,  the  solution  process  will  be  less  likely  to  halve  any  particular  solu¬ 
tion  step  size;  smaller  numbers  would  increase  the  likelihood  of  reduced 
step  sizes.  If  larger  numbers  are  used  for  the  If  and  12  entries,  the  solution 
process  is  more  likely  to  increase  any  particular  solution  step  size;  smaller 
numbers  make  increased  step  sizes  less  likely.  The  user  should  realize  that 
any  increase  in  solution  speed  that  may  result  from  adjustment  of  these  num¬ 
bers  must  necessarily  come  at  the  expense  of  integration  accuracy. 

2.8.5  IMPLICIT  METHOD 

The  integration  methods  described  in  subsections  2. 8. 1,  2. 8. 2  and 
2. 8.  3  are  all  explicit  in  form  and  can  be  used  interchangeably  within  the 
mathematical  formulation  of  SCEPTRE  without  difficulty.  If,  however,  any 
implicit  method  is  to  be  used,  whether  single  step  or  multistep,  an  addi¬ 
tional  computational  step  is  necessary. 
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2. 8.5. 1  Basic  Implicit  Format 


A  generalized  form  of  all  implicit  integration  methods  can  be  written 
as 

P  P 

Y(n+1)  =  J  at  y(n_1)  +  h  2  bi  (74) 

i=0  i=-l 

where  Y  is  the  vector  of  system  state  variables,  Y  is  the  vector  of  state 
variable  derivatives,  n  is  the  step  number,  h  is  the  step  size  and  the  a., 
b.  are  suitably  chosen  constants.  Multistep  methods  are  introduced  if  1 
F>  o.  The  simplest  derivative  form  of  equation  (74)  is 

Y(n+1)  =  Y(n)  +  h  Y(n+1) 


which  is  commonly  referred  to  as  the  implicit  or  backward  Euler  technique. 
If  an  mth  order  system  of  differential  equations  is  to  be  solved  by  this  meth¬ 
od,  the  following  generalized  matrix  equation  will  result: 


dY, 


1  -hg-^  (Yl, . . .  Ym,  t) 


bY 


-hg  (Yl, . . .  Ym,  t) 


r  *1 

.  .  -hg  ^-(Yl,...Ym,t) 
m 

• 

• 

# 

rH 

^  •  •  • 
<1 

.  .  1-hg  bYm  (Yl...Ym,t) 
bY 

m 

lr 

AY 

m 

-F1(Yl,..Ym,t) 

-Fm(Yl,..Ym,t) 


(75) 


The  iteration  implied  in  equation  (75)  is  carried  out  to  convergence  at  each 
time  step.  The  k  superscripts  here  indicate  the  kth  approximation  to  the 
final  value  at  convergence  at  each  step,  g  is  a  constant  that  depends  on  the 
order  of  integration,  and  Fj  is  a  function  of  the  ith  differential  equation, 
the  step  size  and  past  values  of  the  mth  state  variable.  Sparse  matrix  tech¬ 
niques  will  be  applied  to  the  operation  implied  by  equation  (75)  when  large 
problems  are  encountered. 
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2, 8.5.2  The  Jacobian 


The  Y.  terms  in  equation  (75)  are  readily  available  from  the  basic 
program  formulation;  but  the  general  partial  derivative  term 

dY. 

W-,  1  -  i  -  m,  1  -  j  ^  m 

?Yi 

is  not.  To  determine  what  is  reat'y  needed,  it  is  desirable  to  frame  the 
entire  derivation  of  the  sunibolic  Jacobian  in  terms  of  the  general  matrix 
equation 

Y  =  AY  +  BU  +  NU  (76) 

Since  the  desired  quantities  are  the  general 

8Y. 

l 

8Y. 

1 

it  can  be  seen  from  partial  differentiation  of  equation  (76)  that  these  are 
contained  in  the  matrix  A.  Hence  the  convenience  of  the  general  notation 
given  in  equation  (76). 

What  now  follows  is  the  construction  in  symbolic  form  of  the  general  ma¬ 
trix  A  in  terms  of  the  mathematical  formulation  that  was  derived  in  sub¬ 
section  2.3.  Begin  with 


C44^4  -  !4  -  Bi4  l\  +  B24  l2  +  B34  X3  +  B84  J8 


(5) 


and 


V3  "  B34  V4  "  B35  V5  ‘  B36  V6  '  B37  E7 
Substitute  equation  (9)  and  (13)  into  equation  (5)  for 


(4) 


CAA  VA  =  B,/  C1t(-B1AVA-B„EJ+  B24T  Mr1  (-B24V4-B27E? 


'44  4  "  14  "IT ~14*4 


"  B25  R55  B35  h  "  B25  R55  B85  V  +  B34  *3  f  B84  J8  (77) 
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Use  equations  (24),  (10)  (23)  and  (15)  into  equation  (4)  for 


L33  *3  +  L36  ^B36  *3  +  B86  J'8^ 


(78) 


”B34V4  '  B35MG 


-1 


'B25  G22  *D24V4  +  B27E7^  +  B35  *3  f  B85  ll8 

rT 

tv  ; 


_B3C  (L63  +  L66  B36  *  *3  "  B36  L66  BCC  J8  "  B37  E7 


Equation  (74)  can  be  manipulated  to  yield 


<C44  +  B14T  C11B14>  *4  ’  -D14T  C11  B17  E7  '  B24T  Mn'1(D24  V4 


"B27E7  '  B25R55B35  *3  "  B25R55B85  J8 


+  B34T  !3  +  B84T  J8 


(77’) 


Equation  (78)  yields 

J33  +  L36  B36 


f  T  T 

8^10  ^  8pjr  1 J  n  5nr  I  Ir 


'36  63  36  66  "36 


(78’) 


"B34V4  "  B35MG  ^  "B25  G22B24V4  _B25  G22B27E7  +  B35  *3  + 


T 

B85  J8 


)  -  B37E7  - 


L36B86  ^8  "  B36L66B86  ^8 
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If  only  the  coefficients  of  V4,  13,  V4  and  I3  (the  state  variables 
and  the  state  variable  derivatives)  are  retained,  equations  (77')  and  (78’) 
become  considerably  simplified  and  can  be  written  in  matrix  notation  as 


M  0 
c 


o  M. 


-1 


<-B24T  V1  B24>  (B34T  -  B24T  MR  ^  B25R55B35T> 
*B35MG  B25  G22  B24  "  B24^  ^'B35MG  B35  ^ 


(79) 


where  Me  =  C44  +  B14  Cll  B14 


and  Ml  L33  +  L36B36  +  B36Lg3  +  B3gL66B36 


The  matrix  appearing  on  the  right  side  of  equation  (79)  is  now  in  a  form 
corresponding  to  the  first  term  on  the  right  side  of  equation  (76).  This 
matrix  is  the  symbolic  form  of  the  general  A  matrix  that  is  needed  to  use 
implicit  integration  (Ref.  8,  9)  with  the  basic  program  formulation. 


Another  method  is  available  to  construct  the  Jacobian  that  is  based 
on  a  numerical  rather  than  a  symbolic  approach.  Each  time  that  the  Ja¬ 
cobian  is  to  be  reevaluated,  m  calls  are  made  to  SIMUL  8  to  compute 

Yi'  =  $  (Yl’  ....  Y'm,  t) 


Approximations  to  the  desired  partials  are  then  obtained  by 


3  Y. 
1 

3y: 

j 


jv  -  *1 

V,  -  Y. 


(80) 
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Section  III 


SYSTEM  OPERATION 


3.1  INTRODUCTION 

The  SCEPTRE  System  is  a  Fortran  IV  computer  program  written  for 
the  IBM  7090/94  Data  Processing  System.  It  consists  of  two  major  phases. 
The  first,  called  the  Program  Generator,  creates  (on  magnetic  tape)  another 
FORTRAN  IV  computer  program  containing  circuit  equations  for  electrical 
networks.  SCEPTRE  does  this  automatically,  using  the  input  data  describing 
the  circuit  to  be  analyzed.  Both  d-c  steady-state  and  transient  equations  for 
non-linear  circuits  can  be  generated.  The  second  phase,  the  Circuit  Solution 
Executive,  computes  the  circuit  response  by  solving  these  equations  generated 
in  the  FORTRAN  IV  program. 

Operation  of  the  SCEPTRE  System  depends  upon  the  monitor  system 
--D3SYS-- which  controls  its  execution  on  the  IBM  7090/94.  Under  control 
of  IBSYS,  SCEPTRE  executes  in  two  phases  as  separate  job  steps  that  are 
loaded  and  executed  sequentially.  Prior  to  loading,  however,  the  monitor 
system  performs  any  required  FORTRAN  compilation.  Programs  to  be  com¬ 
piled  do  not  have  to  reside  on  the  standard  input  tape  upon  which  the  input 
data  are  stored.  Instead,  the  system  can  be  instructed,  via  system  control 
cards,  to  compile  and  load  programs  from  an  alternate  input  tape.  SCEPTRE 
uses  this  feature  of  IBSYS  as  a  means  of  linking  its  two  job  steps.  This 
is  illustrated  in  the  System  Flow  Diagram,  figure  10. 

3.2  PROGRAM  GENERATOR 


The  Program  Generator  is  an  executive  program  that  controls  the  in¬ 
putting  of  circuit  description  data,  generation  of  a  FORTRAN  IV  subprogram 
for  calculating  circuit  response,  generation  of  circuit  parameter  data,  stor¬ 
age  of  circuit  models  on  library  tapes,  restart  of  discontinued  runs,  and  re¬ 
outputting  of  computed  results.  Each  of  these  six  program  tasks  as  well  as 
the  Program  Generator  (EXEC1)  are  described  in  the  flow  diagrams,  figures 
11  through  16. 

3.2.1  CIRCUIT  DESCRIPTION  PROCESSOR 

The  SCEPTRE  circuit  description  language  is  used  to  describe  elec¬ 
trical  networks.  This  application -oriented  language  is  powerful,  easy  to 
learn  and  use,  and  nearly  format  free.  With  it,  circuits  composed  of  funda¬ 
mental  circuit  elements  and  prestored  circuit  models  can  be  described.  The 
types  of  fundamental  electrical  characteristics  allowed  are  resistance. 
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Figure  11.  Program  Generator  (EXEC  1)  Flow  Diagram 
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Figure  12.  Circuit  Description  Processor  (CDPRS)  Flow  Diagram 


WRTTPE _ 

Pr  inf  -out  Final 
Permanent  Library 
If  'Print'  Requested 


V.'RTTPE _ 

Prinl-ouf  Final 
Temporary  Library, 
If  any 


Enter 

MDEDIT 


Read  Permanent 
Library  Index 


Read  a  Cord 
from  Tape 


Mo  den 
Header 
.Card  ?v 


Get  Model  Name 


Permanent 
Library 
Model  ? 


Add  Model  to 
Permanent  Llbrory 
Index 


Add  Data  Cord 
To  Indicated 
Model  Library 


to 

Temporary 
Library  Index 


Figure  13.  Model  Editor  (MDEDIT)  Flow  Diagram 
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Figure  15.  Rerun  Processor  and  Data  Generator  (RRPRO)  Flow  Diagran 
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capacitance,  inductance,  current  and  voltage  sources,  and  mutual  inductance. 
Using  the  same  components  and  circuit  description  langauge,  equivalent  cir¬ 
cuit  models  of  devices  such  as  transistors  and  diodes  can  be  described  and 
stored  on  magnetic  tape  for  future  use.  When  a  stored  model  is  referenced 
in  a  circuit  description,  it  is  located  on  the  specified  library  tape  and  its 
elements  appropriately  substituted  into  the  circuit  being  described. 

Often,  when  a  model  is  called  out  in  describing  a  circuit,  the  desired 
parameter  values  are  different  from  those  originally  specified.  Rather  than 
permanently  change  the  model  stored  on  tape,  the  SCEPTRE  circuit  descrip¬ 
tion  language  allows  n  odel  parameter  valups  to  be  changed  when  model  ele¬ 
ment  substitution  takes  place. 

3.2.2  MODEL  EDITOR 

Using  the  SCEPTRE  circuit  description  language,  circuit  models  con¬ 
structed  by  the  user  can  be  stored  on  a  model  library  tape  by  the  Model 
Editor  program.  One  permanent  library  tape  and  one  temporary  library  tape 
are  provided  for  storing  models.  Circuit  models  may  be  any  arbitrary  n- 
terminal  configuration  of  the  allowed  fundamental  circuit  elements.  Any 
model  so  described  can  be  stored  on  either  library  tape  using  the  Model  Edi¬ 
tor  program. 

3.2.3  CIRCUIT  EQUATION  GENERATOR 

After  the  description  of  the  circuit  to  be  analyzed  has  been  reconstructed 
in  memory  by  the  Circuit  Description  Processor,  the  FORTRAN  IV  sub¬ 
program  called  SIMUL8  is  created.  It  contains  either  d-c  steady-state  equa¬ 
tions,  transient  solution  equations,  or  both,  depending  on  the  type  of  analysis 
requested  by  the  user.  These  equations  are  based  on  the  formulation  pre¬ 
sented  inSection  II  of  this  report.  The  SIMUL8  program  is  written  on  mag¬ 
netic  tape  (PROGRAM  SAVE  TAPE)  and  stored  until  the  second  phase  of 
SCEPTRE  operation,  when  it  will  be  compiled  and  executed.  The  manner  in 
which  the  circuit  equations  are  written  provides  a  very  efficient  computation 
of  the  circuit  response  during  execution,  thus  allowing  the  response  of  a 
large  circuit  to  be  calculated  rapidly. 

3.2.4  DATA  GENERATOR  AND  RERUN  PROCESSOR 

Essentially  the  SIMUL8  program  contains  only  the  circuit  equations  for 
the  network  under  investigation.  The  parameter  or  component  values  of  the 
network  are  stored  separately  as  input  data  to  the  SIMUL8  program.  The 
Data  Generator  program  organizes  and  stores  the  circuit  data  on  the  PRO¬ 
GRAM  SAVE  TAPE. 
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Once  a  circuit  has  been  described  to  SCEPTRE,  multiple  or  repeated 
circuit  solutions  can  be  run  by  changing  parameter  values  between  runs. 

The  circuit  description  language  is  used  to  specify  the  number  of  repeated 
runs,  each  known  as  a  rerun,  and  the  changes  in  parameter  values  desired 
for  each  rerun.  The  Rerun  Processor  interprets  and  processes  this  informa¬ 
tion.  The  Da^a  Generator  then  creates  and  stores  one  block  of  circuit  data 
for  each  rerun  on  the  PROGRAM  SAVE  TAPE.  Since  only  parameter  values 
are  changed  between  reruns  (i.  e. .  no  topological  changes),  the  SIMUL8  pro¬ 
gram  containing  the  original  circuit  equations  is  simply  re-executed  during 
the  solution  phase  using  the  data  created  for  each  rerun. 

3.2.5  CONTINUE  PROCESSOR 

The  Continue  Processor  allows  previously  terminated  transient  solution 
runs  to  be  restarted  and  continued.  In  addition,  it  is  necessary  that  the 
PROGRAM  SAVE  TAPE  from  the  original  transient  solution  run  be  available. 
This  tape  contains,  in  addition  to  the  SIMUL8  program,  all  of  the  parameter 
values  and  data  required  to  continue  the  transient  solution.  The  Continue 
Processor  also  allows  certain  run  control  parameters  to  be  changed  and  used 
for  the  continuation  of  the  original  transient  run. 

3.2.6  RE-OUTPUT  PROCESSOR 

The  SCEPTRE  user  may  wish  to  re-create  both  the  printed  and  plotted 
outputs  for  a  particular  transient  analysis  run.  This  can  be  done,  providing 
the  OUTPUT  SAVE  TAPE  is  saved  from  the  original  solution  run. by  using  the 
Re-Output  Processor.  Several  changes  in  the  output  can  be  made  during  a 
Re-Output  run.  For  example: 

•  Ou  tput  labels  can  be  changed 

•  New  quantities  can  be  plotted 

•  The  order  in  which  output  quantities  are  printed  out  and  plotted 
can  be  changed 

•  Printing  and  plotting  of  output  quantities  can  be  suppressed. 

3.3  SOLUTION  EXECUTIVE 


In  the  second  phase  of  SCEPTRE  operation,  the  FORTRAN  IV  subpro¬ 
gram  SIMUL8  is  compiled  and  executed.  During  the  loading  of  this  job  step, 
the  IBSYS  system  is  instructed  to  read  the  SEMUL8  program  from  the  PRO¬ 
GRAM  SAVE  TAPE,  compile  it,  and  then  load  it  along  with  the  other  pro¬ 
grams  being  loaded. 
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Execution  commences,  as  shown  in  figure  16,  with  the  reading  of  the 
run  control  parameters  from  the  PROGRAM  SAVE  TAPE.  Then  the  SIMUL8 
program  is  called.  It  reads  the  remaining  circuit  parameter  values  stored 
on  die  PROGRAM  SAVE  TAPE  and  then  enters  the  circuit  solution  equations. 
The  circuit  equations  are  solved,  thus  computing  the  state-variable  deriva- 
t  ves.  Calling  the  selected  numerical  integration  routine  then  produces  new 
values  of  the  state  variables  for  the  next  time  step.  At  the  conclusion  of  each 
successful  integration  step,  the  requested  output  quantities  are  buffered  in 
memory.  When  the  buffer  is  full,  it  is  written  on  the  OUTPUT  SAVE  TAPE 
in  FORTRAN  IV  binary  format.  After  the  circuit  solution  is  complete,  con¬ 
trol  is  returned  to  the  Solution  Executive  Program,  whereupon  the  contents 
of  the  OUTPUT  SAVE  TAPE  (computed  results)  are  re -formatted  in  lists  and 
graphs  and  stored  on  the  SYSTEM  OUTPUT  TAPE  for  peripheral  processing. 

If  any  reruns  were  requested,  control  is  ti->c-n  returned  to  SIMUL8  for 
re-execution  of  the  solution  phase  and  subsequent  outputting.  This  re¬ 
cycling  is  continued  until  all  circuit  reruns  have  been  processed. 

At  the  conclusion  of  each  transient  solution,  the  critical  parameter 
values  and  data  are  stored  on  the  PROGRAM  SAVE  TAPE.  Thus,  if  the  tape 
is  saved  at  the  end  of  the  run,  the  soluti  a  can  be  continued  at  some  future 
time.  Use  of  the  Continue  Processor  allows  previously  terminated  solution 
runs  to  be  restarted  and  continued  with  changes  in  the  run  control  parameters. 

Transient  solution  runs  may  be  terminated  or  simply  saved  by  the  com¬ 
puter  operator  by  depressing  sense  switch  No.  6  on  the  7090/94  console. 

Using  this  feature,  a  PROGRAM  SAVE  TAPE  could,  for  example,  be  gener¬ 
ated  every  15  minutes  on  long  running  transient  solutions,  eliminating  the 
need  for  repeating  previous  calculations  in  the  event  of  an  abnormal  run 
termination. 

If,  at  the  conclusion  of  a  transient  analysis  run,  the  OUTPUT  SAVE 
TAPE  is  saved,  the  Re-Output  Processor  may  be  used  to  reproduce  lists  and 
graphs  of  any  of  the  circuit  quantities  originally  requested  for  output.  In 
addition,  graphs  of  variables  plotted  against  variables  other  than  time  which 
may  not  have  been  requested  originally  can  be  conveniently  produced. 
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Appendix  I 

B  MATRIX  DERIVATION 


The  derh'ation  of  the  general  B  matrix  that  expresses  link  voltages  in 
terms  of  tree  branch  voltages  and  tree  branch  currents  in  terms  of  link 
currents  is  <ts  follows; 


Two  fundamental  incidence  matrices  that  arise  from  network  topology 
theory  will  be  called  the  Q  and  T  matrices  here.  The  fundamental  cut  set 
matrix,  Q  =  [q^  is  a  matrix  containing  (n-1)  rows  and  b  columns  for  a  net¬ 
work  containing  n  nodes  and  b  elements,  where; 


+1  if  the  i**1  fundamental  cut  set  direction  coincides  with  the  refer¬ 
ence  direction  of  the  j*h  element 


-1  if  the  i**1  fundamental  cut  set  direction  is  in  opposition  to  the 
reference  direction  of  the  ji*1  element 


0  if  the  i*-*1  fundamental  cut  set  does  not  include  the  j111  element 


If  the  elements  are  properly  ordered,  it  is  always  true  that 

Q  =  [-BT  u] 


where  the  columns  of  the  unit  matrix  U  correspond  to  the  tree  branch  ele 
ments. 


The  fundamental  circuit  matrix,  T  =  C  t  j j  J  is  a  matrix  containing  m 
rows  and  b  columns  for  a  network  containing  m  independent  loops  and  b  ele¬ 
ments,  where: 


t..  -  fl  if  the  i**1  independent  loop  direction  coincides  with  the  refer- 
1J  once  direction  of  the  ji*1  element 


ty  -  -1  if  the  i^h  independent  loop  direction  opposes  the'  refer¬ 
ence  direction  of  the  j**1  element 

t . -  --  0  if  the  i*h  independent  loop  does  not  include  the  j1*1  element 
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If  the  elements  are  properly  ordered,  it  is  always  true  that 
T  =  [U  B] 

where  the  columns  of  the  unit  matrix  U  correspond  to  the  network  links. 
Since  it  is  always  true  that 
Q  I  b  =  0.  T  Vb  -  0 
direct  substitution  yields 


-bt  u 

jL  ' 

r  n 

=  0  and  UBl 

vL  ‘ 

L  J 

L  ‘tbJ 

L  J 

_  VTB. 

Expansion  of  these  relations  gets 

ITB  ^  B  !L  a^d  VL  77  -B  Vxn 

so  that  the  tree  branch  currents  may  be  expressed  in  terms  of  the  link  cur 
rents  and  the  link  voltages  may  be  expressed  in  terms  of  the  tree  branch 
voltages  through  the  B  matrix. 
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Appendix  TI 

DIODE  REPRESENTATION  IN  THE  INITIAL-CONDITIONS  PROGRAM 


The  current  in  a  diode  or  transistor  junction  at  a  point  on  its  voltage  vs 
current  (V/I)  characteristic  may  be  separated  into  two  components  as 
J  -  GjVj+Q.  Consider  the  typical  diode  curve  below. 


Angle  aj  =  angle  03  is  enclosed  by  the  slope  of  the  diode  characteristic 
at  any  point  and  the  horizontal  at  that  point.  E'  is  an  offset  voltage  that  marks 
the  intersection  of  a  continuation  of  the  slope  Gj  and  the  line  J  =  0. 

From  the  figure: 

tan  a  x  =  tan  a  2  =  Gj  =  yJT'g7 

or  J  :■  G  j  Vj  -  Gj  E' 

or  J  G  j  V  j  t  Q  if  Q  is  defined  as  -  G  jE' 
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Appendix  III 

BASIC  NEWTON -RAPHSON  METHOD 


Given  a  single  algebraic  or  transcendental  equation  of  the  form  F(x)  -  0, 
that  is  single  valued  and  differentiable  in  the  domain  of  interest,  a  Newton - 
Raphson  procedure  can  be  constructed  using 


F(x)  +  — ^  Ax  =  0 

d  X 

An  initial  x  -  xq  may  be  assumed  and  F(xq), 


9  F(xq) 
8  x 


determined. 


Then  Ax  = 


-F(xo) 

d  X 


and  x 


1 


'0 


Ax 


The  procedure  is  repetitive  until 


F(xn+i)  -  Ft*,,) 


<  z 


where  z  is  some  specified  convergence  criteria.  When  this  last  relation  is 
satisfied,  the  process  is  said  to  have  converged.  Extension  to  systems  of 
equations  adds  no  complications. 


Appendix  IV 

EQUIVALENCE  OF  EQUATIONS  (44)  AND  (44x) 


To  show  the  equivalence  of  equations  (44)  and  x  (44 f).  it  will  be 
sufficient  to  show  that 


>9' 

-1 

'Fl  (V9,  V2,  V5)‘ 

-1 

-  B97  E7 

V2 

-  N 

f2  (v9i  v2.  V5) 

-  1Z] 

-1327  F7 

V5 

F3  (Vg,  V2,  V5) 

[b95T+  b05  ^  |  *  B85  r,I8 

or 


’v9' 

r 

-D97  E7 

F,  (Vg.  V2,  V5) 

•x 

v2 

=  [z]  < 

-B27  E7 

F2  (Vg,  V2,  V5) 

v5 

L  - 

|b95T+  B05Tq1  ^9  +  bS5  J8 

f3  <v9’  v2-  V5> 

J 

or 


V9" 

-B97  E7 

Fl  (Vg.  V2.  V5) 

V2 

= 

-B27  E7 

F2  (Vg.  V2.  Vg) 

V5 

|b95T+  B05Ta]  Qg  +  B851J8 

F3  (Vg,  V2.  V5) 

The  left  side  of  the  alx»ve  equation  expands  into 


V9  *  D95  V5 
V2  *  B25  V5 


1  I 

B95  4  B0S  a  G99V9  -  B25TG22V2  +  G55V5 


59 


4 


and  the  right  side,  upon  .ub.il, r„oni„ms  (37,  th,,,^  0#)  llel,!nl(.s 


V9  *  B95  '  h 

V2  ’  ^  V, 


'  [b95  +  b05^j]  G99  V9  -  C.22  V'2  -  G53  V.s 

which  shows  the  equivalence 
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